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Abstract The classical macroscopic chemotaxis equations have previously been 
derived from an individual-based description of the tactic response of cells that use 
a "run-and-tumble" strategy in response to environmental cues 1 17 18 1. Here we 
derive macroscopic equations for the more complex type of behavioral response 
characteristic of crawling cells, which detect a signal, extract directional informa- 
tion from a scalar concentration field, and change their motile behavior accord- 
ingly. We present several models of increasing complexity for which the deriva- 
tion of population-level equations is possible, and we show how experimentally- 
measured statistics can be obtained from the transport equation formalism. We also 
show that amoeboid cells that do not adapt to constant signals can still aggregate 
in steady gradients, but not in response to periodic waves. This is in contrast to the 
case of cells that use a "run-and-tumble" strategy, where adaptation is essential. 

1. Introduction 

Motile organisms sense their environment and can respond to it by (i) directed 
movement toward or away from a signal, which is called taxis, (ii) by changing 
their speed of movement and/or frequency of turning, which is called kinesis, or 
(iii) by a combination of these. Usually these responses are both called taxes, and 
we adopt this convention here. Taxis involves three major components: (i) an ex- 
ternal signal, (ii) signal transduction machinery for transducing the external signal 
into an internal signal, and (iii) internal components that respond to the trans- 
duced signal and lead to changes in the pattern of motility. In order to move away 
from noxious substances (repellents) or toward food sources (attractants) organ- 
isms must extract directional information from an extracellular scalar field, and 
there are two distinct strategies that are used to do this. A simple paradigm will 
illustrate these. 
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Suppose that one is close enough to a bakery to detect the odors, but cannot 
see the bakery. To find it, one strategy is to use sensors at the end of each arm that 
measure the difference in the signal at the current location and use the difference to 
decide on a direction. Clearly humans do not use this strategy, but instead, execute 
the "bakery walk", which is to take a sniff and judge the signal intensity at the 
present location, take a step and another sniff, compare the signals, and from the 
comparison decide on the next step. 

The first strategy is used by amoeboid cells (cells which move by crawling 
through their environment), which have receptors on the cell membrane and are 
large enough to detect typical differences in the signal over their body length. 
Small cells such as bacteria cannot effectively make a "two-point in space" mea- 
surement over their body length, and therefore they adopt the second strategy and 
measure the temporal variation in the signal as they move through the external 
field. In either case, an important consideration in understanding population-level 
behavior is whether or not the individual merely detects the signal and responds 
to it, or whether the individual alters it as well, for example by consuming it or 
by amplifying it so as to relay the signal. In the former case there is no feedback 
from the local density of individuals to the external field, but when the individ- 
ual produces or degrades the signal, there is coupling between the local density of 
individuals and the intensity of the signal. The latter occurs, for example, when 
individuals move toward a signal from neighboring cells and relay the signal as 
well, as in the aggregation of the cellular slime mold Dictyostelium discoideum 
(Dd). 

One of the best-characterized systems that adopts the "bakery walk" strategy 
is the flagellated bacterium E. coli, for which the signal transduction machinery 
is well characterized |5|. E. coli alternates between a more or less linear motion 
called a run and a highly erratic motion called tumbling, which produces little 
translocation but reorients the cell. Run times are typically much longer than the 
tumbling time, and when bacteria move in a favorable direction (i.e., either in 
the direction of foodstuffs or away from noxious substances), the run times are 
increased further. Since these bacteria are too small to detect spatial differences 
in the concentration of an attractant on the scale of a cell length, they choose a 
new direction essentially at random at the end of a tumble, although it has some 
bias in the direction of the preceding run |4|. The effect of alternating these two 
modes of behavior and, in particular, of increasing the run length when moving in a 
favorable direction, is that a bacterium executes a three-dimensional random walk 
with drift in the favorable direction, when observed on a sufficiently long time 
scale 1 3 30 1. In addition, these bacteria adapt to constant signal levels and in effect 
only alter the run length in response to changes in extracellular signals. Models for 
signal transduction and adaptation in this system has been developed |46 2\, and a 
simplified version of the first model has been incorporated into a population-level 
description of behavior 11711 81 . The latter analysis shows how parameters that 
characterize signal transduction and response in individual cells are embedded in 
the macroscopic sensitivity x m the macroscopic chemotaxis equation described 
later. Having the bacterial example in mind, we will call the "bakery walk" strategy 
as a "run-and-tumble" strategy in what follows. 
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The directed motion of amoeboid cells (e.g. Dd or leukocytes), which is cru- 
cial in embryonic development, wound repair, the immune response to bacterial 
invasion, and tumor formation and metastasis, is much more complicated than 
bacterial motion. Cells detect extracellular chemical and mechanical signals via 
membrane receptors, and these trigger signal transduction cascades that produce 
intracellular signals. Small differences in the extracellular signal over the cell are 
amplified into large end-to-end intracellular differences that control the motile ma- 
chinery of the cell and thereby determine the spatial localization of contact sites 
with the substrate and the sites of force-generation needed to produce directed 
motion 1 40 7 1. Movement of Dd and other amoeboid cells involves at least four 
different stages 13 61431 . (1) Cells first extend localized protrusions at the leading 
edge, which take the form of lamellipodia, filopodia or pseudopodia. (2) Not all 
protrusions are persistent, in that they must anchor to the substrate or to another 
cell in order for the remainder of the cell to follow 1 44 1 . Protrusions are stabilized 
by formation of adhesive complexes, which serve as sites for molecular signaling 
and also transmit mechanical force to the substrate. (3) Next, in fibroblasts acto- 
myosin filaments contract at the front of the cell and pull the cell body toward the 
protrusion, whereas in Dd, contraction is at the rear and the cytoplasm is squeezed 
forward. (4) Finally cells detach the adhesive contacts at the rear, allowing the tail 
of the cell to follow the main cell body. In Dd the adhesive contacts are relatively 
weak and the cells move rapidly (~ 20/xm/min), whereas in fibroblasts they are 
very strong and cells move slowly. The coordination and control of this complex 
process of direction sensing, amplification of spatial differences in the signal, as- 
sembly of the motile machinery, and control of the attachment to the substratum 
involves numerous molecules whose spatial distribution serves to distinguish the 
front from the rear of the cell, and whose temporal expression is tightly controlled. 
In addition, Dd cells adapt to the mean extracellular signal level |40|. 

Dd is a widely-used model system for studying signal transduction, chemo- 
taxis, and cell motility. Dd uses cAMP as a messenger for signaling initiated by 
pacemaker cells to control cell movement in various stages of development (re- 
viewed in 1 39 1). In the absence of cAMP stimuli Dd cells extend pseudopods 
in more-or-less random directions, although not strictly so since formation of a 
pseudopod inhibits formation of another one nearby for some time. Aggregation- 
competent cells respond to cAMP stimuli by suppressing existing pseudopods and 
rounding up (the "cringe response"), which occurs within about 20 sees and lasts 
about 30 sees |8 1. Under uniform elevation of the ambient cAMP this is followed 
by extension of pseudopods in various directions and an increase in the motility 
1531541 . However, one pseudopod usually dominates, even under uniform stimu- 
lation. A localized application of cAMP elicits the "cringe response" followed by 
a localized extension of a pseudopod near the point of application of the stimulus 
1471 . How the cell determines the direction in which the signal is largest, and how 
it organizes the motile machinery to polarize and move in that direction, are ma- 
jor questions from both the experimental and theoretical viewpoint. Since cAMP 
receptors remain uniformly distributed around the cell membrane during a tac- 
tic response, receptor localization or aggregation is not part of the response |27|. 
Well-polarized cells are able to detect and respond to chemoattractant gradients 
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with as little as a 2% concentration difference between the anterior and posterior 
of the cell 1 34 1 . Directional changes of a shallow gradient induce polarized cells to 
turn on a time scale of 2-3 seconds |23 1, whereas large changes lead to large-scale 
disassembly of motile components and creation of a new "leading edge" directed 
toward the stimulus 1 22 1. Polarity is labile in cells starved for short periods in that 
cells can rapidly change their leading edge when the stimulus is moved |47|. 

There are a number of models for how cells extract directional information 
from the cAMP field. Fisher et al. 1198 suggest that directional information is 
obtained by the extension of pseudopods bearing cAMP receptors and that sensing 
the temporal change experienced by a receptor is equivalent to sensing the spatial 
gradient. However, Dd cells contain a cAMP-degrading enzyme on their surface, 
and it has been shown that as a result, the cAMP concentration increases in all 
directions normal to the cell surface 1 10 1. Furthermore, more recent experiments 
show that cells in a steady gradient can polarize in the direction of the gradient 
without extending pseudopods |40|. Thus cells must rely entirely on differences 
in the signal across the cell body for orientation. Moreover, the timing between 
different components of the response is critical, because a cell must decide how to 
move before it begins to relay the signal. Analysis of a model for the cAMP relay 
pathway shows that a cell experiences a significant difference in the front-to-back 
ratio of cAMP when a neighboring cell begins to signal 1 10 1, which demonstrates 
that sufficient end-to-end differences for reliable orientation can be generated for 
typical extracellular signals. An activator-inhibitor model for an amplification step 
in chemotactically sensitive cells was also postulated 1351 . Amplification of small 
external differences involves a Turing instability in the activator-inhibitor system, 
coupled to a slower inactivator that suppresses the primary activation. While this 
model reproduces some of the observed behavior, there is no biochemical basis for 
it; it is purely hypothetical and omits some of the major known processes. A model 
that takes into account some of the known biochemical steps has been proposed 
more recently ITTI . 

The objective of this paper is to derive equations for the population-level be- 
havior of amoeboid cells such as Dd or leukocytes that incorporate details about the 
individual-based response to signals. We present several models with the increased 
complexity for which the derivation of population-level equations is possible. We 
show how experimentally-measured statistics can be obtained from the transport 
equation formalism. The paper is organized as follows. We discuss the classical 
chemotaxis description and summarize the state of the art of the derivation of 
macroscopic equations and population-level statistics from individual-based mod- 
els in the remainder of this section. In Section|2] we establish the general setup for 
models of amoeboid cells and we present individual-based models which capture 
the essential behavioral responses of eukaryotic cells. In Section[5]we derive the 
macroscopic moment equations from the microscopic model and the dependence 
of the mean speed on the signal strength is studied. Finally, we provide conclusions 
and the discussion of the presented approaches in Section[4] 
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1.1. Macroscopic descriptions of chemotaxis 

The simplest description of cell movement in the presence of both diffusive and 
tactic components results by postulating that the flux of cells j is given by 

j = -DVn + nu c , (1) 

where n is the density of cells, u c is the macroscopic chemotactic velocity and D 
is the difusion constant. The taxis is positive or negative according as u c is parallel 
or anti-parallel to the direction of increase of the chemotactic substance S. Keller 
and Segel [28 1 postulated that the chemotactic velocity is given by u c = x(S)VS 
and then can be written as 

j = -DVn + n X (S)VS (2) 

where x(S) is called the chemotactic sensitivity. In the absence of cell division or 
death the resulting conservation equation for the cell density n(x, t) is 

3 TL 

— = V • (DVn - nx(S)VS) (3) 

and this is called a classical chemotaxis equation. Unless the distribution of the 
chemotactic substance is fixed, l|3} is coupled to an evolution equation for this 
substance, and perhaps other governing variables. 

Other phenomenological approaches to the derivation of the chemotactic ve- 
locity have been taken. For example, by approaching taxis from a mechanical point 
of view, Pate and Othmer 1 41 1 derived the velocity in terms of forces exerted by an 
amoeboid cell. Starting from Newton's law, neglecting inertial effects, and assum- 
ing that the motive force exerted by a cell is a function of the attractant concentra- 
tion, they showed how the chemotactic sensitivity is related to the rate of change of 
the force with attractant concentration. In this formulation the dependence of the 
flux on the gradient of the attractant arises from the difference in the force exerted 
in different directions due to different attractant concentrations. Experimental sup- 
port for this comes from work of |51 1, who show that as many pseudopods are 
produced down-gradient as up, but those up-gradient are more successful in gen- 
erating cell movement. We shall use a version of the mechanical approach to taxis 
in a model described in the following section. 

The first derivation that directly relates the chemotactic velocity to properties 
of individual cells is due to Patlak |42|, who used kinetic theory arguments to 
express u c in terms of averages of the velocities and run times of individual cells. 
This approach was extended by Alt 1 1 1, who showed that for a class of receptor- 
based models the flux is approximately given by (0. These approaches are based 
on velocity-jump processes, which lead to transport equations of the form 

—p(xL,v,t)+v-Vp(x,v,t)=-\pCx,v,t)-\-\ T(v,v')p(x,v',t)dv'. (4) 
at Jv 

where p(x, v,i) is the density of cells at position x 6 ft C R n , moving with 
velocity v £ V C R" at time t > 0, A is the turning rate and kernel T(v, v') 
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gives the probability of a change in velocity from v' to v, given that a reorienta- 
tion occurs 1 37 1. External signals enter either through a direct effect on the turning 
rate A and the turning kernel T, or indirectly via internal variables that reflect the 
external signal and in turn influence A and/or T. The first case arises when exper- 
imental results are used to directly estimate parameters in the equation |20|, but 
the latter approach is more fundamental. The reduction of ® to the macroscopic 
chemotaxis equations for the first case is done in |24 38 1 and |6|. 

Some statistics of the density distribution in the first case, wherein the external 
field modifies the turning kernel or turning rate directly, can easily be derived and 
used to interpret experimental data. To outline the procedure, we consider two- 
dimensional motion of amoeboid cells in a constant chemotactic gradient directed 
along the positive x% axis of the plane, i.e. 

VS = ||V5|| ei, where we denoted ei = [1, 0]. (5) 

Moreover, we assume that the gradient only influences the turn angle distribu- 
tion T; details of the procedure are given in 1371 . We assume for simplicity that 
the individuals move with a constant speed s, i.e. a velocity of an individual can 
be expressed as v(0) = s[cos(0), sin(</>)] where <f> E [0, 2ii). We assume that 
T(v, v') = T(<fi, (/>') is the sum of a symmetric probability distribution h(<f>) and 
a bias term k{<f>) that results from the gradient of the chemotactic substance. Since 
the gradient is directed along the positive X\ axis, we assume that the bias is sym- 
metric about (p = and takes its maximum there. Thus we write T(cf>, (/>') — 
h{(f) — 4>') + k((f>) where h and k are normalized as follows. 

27T r2TT 

h(<f>)d<j> = 1 / fe(0)d0 = (6) 

Jo 

Let p(x, <fi, t) be the density of cells at position x E M 2 , moving with velocity 
v(cf>) = s[cos(4>), sm(<p)], <j> E [0, 2tt), at time t > 0. The statistics of interest are 
the mean location of cells X(t), their mean squared displacement V 2 (t), and their 
mean velocity V(£), which are defined as follows. 

Jr 2 Jo 

^0 JR2 Jo 



Jk 2 Jo 



(0))p(x,0,i) d^dx, 



where No is the total number of individuals present and B[t) is an auxiliary vari- 
able that is needed in the analysis. Two further quantities that arise naturally are 
the taxis coefficient x, which is analogous to the chemotactic sensitivity defined 
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earlier because it measures the response to a directional signal, and the persistence 
index tpd- These are defined as 

X= cos 4>d(j) and -0d = 2 / h(4>) cos 4>&4>- (7) 

Jo Jo 

The persistence index measures the tendency of a cell to continue in the current 
direction. Since we have assumed that the speed is constant, we must also assume 
that x ana 4>d satisfy the relation % < 1 — tpd, for otherwise the former assumption 
is violated (cf. JTol). 

One can now show, by taking moments of 0, using and symmetries of h 
and k, that the moments satisfy the following evolution equations [37 1. 

f=V ^ = -AoV + Ax, ei (8) 

= 2B ? = - *oB + \ X sX 1 (9) 

at at 

where Ao = A(l — tpd)- The solution of l|8} subject to zero initial data is 

X(t) = sC I (t--j-^l-e- x ° t )^e 1 , V(t) = sC I (l-e- x ° t )e 1 (10) 

where Ci = x/(l — V'd) is sometimes called the chemotropism index. Thus the 
mean velocity of cell movement is parallel to the direction of the chemotactic 
gradient and approaches Voo = s CiB\ as t — > oo. Thus the asymptotic mean 
speed is the cell speed decreased by the factor Cj. 

A measure of the fluctuations of the cell path around the expected value is 
provided by the mean square deviation, which is defined as 

a2{ ^ = Trl 1^ l|x-X(t)|| 2 p(x, ( /»,i)d0dx = P 2 (i)- ||X(t)|| 2 . (11) 

Using (|8j - (|9}, one also finds a differential equation for a 2 . Solving this equation, 
we find 

„>.,£{(, («<-_,)} as 

and from this one can extract the diffusion coefficient as 

2s 2 

Therefore if the effect of an external gradient can be quantified experimentally and 
represented as the distribution k, the macroscopic diffusion coefficient, the persis- 
tence index, and the chemotactic sensitivity can be computed from measurements 
of the mean displacement, the asymptotic speed and the mean-squared displace- 
ment. 

However, it is not as straightforward to derive directly the macroscopic evo- 
lution equations based on detailed models of signal transduction and response. 
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Suppose that the internal dynamics that describe signal detection, transduction, 
processing and response are described by the system 

^ = f(y,S) (12) 

where y € M m is the vector of internal variables and S is the chemotactic sub- 
stance (S is extracellular cAMP for Dd aggregation). Models that describe the 
cAMP transduction pathway exist [33 48 49 1, but for describing chemotaxis one 
would have to formulate a more detailed model. The form of this system can be 
very general but it should always have the "adaptive" property that the steady-state 
value (corresponding to the constant stimulus) of the appropriate internal variable 
(the "response regulator") is independent of the absolute value of the stimulus, and 
that the steady state is globally attracting with respect to the positive cone of K m . 

We showed earlier that for non-interacting walkers the internal dynamics can 
be incorporated in the transport equation as follows 1171 . Let p(x, v, y, t) be the 
density of individuals in a (2N + m)— dimensional phase space with coordinates 
[x, v, y], where x e R N is the position of a cell, v e V C l w is its velocity and 
y G Y C M m is its internal state, which evolves according to (I12t . The evolution 
of p is governed by the transport equation 

^ + V x • wp + V y • ip = -\{y)p + j A(y)T(v, v', y)p(x, v',y, <)dv' (13) 

where, as before, we assume that the random velocity changes are the result of a 
Poisson process of intensity X(y). The kernel T(v, v' y) gives the probability of 
a change in velocity from v' to v, given that a reorientation occurs. The kernel T 
is non-negative and satisfies the normalization condition f v T(v, v',y)dv = 1. 
To connect this with the chemotaxis equation (0, we have to derive an evolution 
equation for the macroscopic density of individuals 

n(x,t)= / / p(x, v, y, t)dvdy. (14) 

J Y JV 

The problem turns out to be tractable for systems that execute "run-and-tumble" 
motion, such as E. coli. To illustrate this, assume for simplicity that the motion is 
restricted to ID, the signal is time-independent, the speed s is constant, and the 
turning phase is neglected; the general cases are treated elsewhere I17I18I . Letp + 
(resp. p~) be the density of individuals moving to the right (resp. left). Then Jl 31 
leads to a telegraph process described by the hyperbolic system 



op Op y—^ o 



s- 



dt dx ^ dy. 



E t~ [ /i(y ' s)p ^ = A(y) [ p+ - p ^ ■ (16) 
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The essential components of the internal dynamics in the bacterial context are fast 
excitation, followed by slower adaptation and return to the basal turning rate, and 
these aspects are captured in the system [ 39 1 



dyi = 9(S(x)) - (yi+y 2 ) and dm _ g(S(x)) - y 2 
dt r e dt r a 

Here g encodes the first step of signal transduction, S is the chemoattractant, and 
r e and r a are time constants for excitation and adaptation, respectively. The com- 
ponent yi adapts perfectly to constant stimuli, i.e., the steady state response is 
independent of the magnitude of the stimulus S. To obtain a macroscopic limit 
equation for the total density n(x, t) we incorporate the variables yi into the state 
and derive a system of four moment equations for various densities and fluxes 11171 . 
Assuming that the turning rate has the form X(y) = Xq — by%, for Ao > 0, b > 0, 
we show that this system reduces to the classical chemotaxis equation for large 
times 



On d 
dt dx V 2Aq dx 



bs 2 T a g'(S(x)) 



A (l + 2Aor Q )(l + 2Aor e ) 



S'(x)nj (18) 



where the chemotactic sensitivity is given explicitly in terms of parameters that 
characterize signal transduction and response. We have only used the simplified 
dynamics (I17> to obtain the macroscopic chemotactic sensitivity, but this model 
captures the essential aspects for bacterial taxis 1461 171 . An open problem is how 
one extracts the elementary processes of excitation and adaptation from a complex 
network of the type used for signal transduction in E. coli. Finally, let us note that 
the global existence results for Jl 3i which is coupled with the evolution equation 
for the extracellular signal were recently given in 1141 . 

Equation Jl 81 was derived for cells such as bacteria, that use the "run-and- 
tumble" strategy, and our objective in this paper is to attempt a similar reduction 
of the transport equation to a chemotaxis equation for more complex amoeboid 
eukaryotic cells. In the following section we introduce the general setup for study- 
ing amoeboid taxis. Then we study several "caricature" or "cartoon" models for 
amoeboid chemotaxis with the objective of deriving macroscopic population-level 
equations in each case. We start with a model which can capture interesting fea- 
tures of eukaryotic motility without introducing additional internal state variables, 
and then add internal state variables to the model. 



2. Amoeboid taxis with internal variables 

A fundamental assumption in the use of velocity-jump processes 1 37 1 to describe 
cell motion is that the jumps are instantaneous, and therefore the forces are Dirac 
distributions. This approximates the case in which very large forces act over very 
short time intervals, and even if one incorporates a resting or tumbling phase, as 
was done in |38|, the macroscopic description of motion is unchanged. This is 
appropriate for the analysis of bacterial motion (and other systems that use a "run- 
and-tumble" strategy), as summarized above, since the effect of the external signal 
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is to change the rotational behavior of the fiagella, and not, so far as it is under- 
stood, to affect the force generation mechanism itself. However, the situation is 
very different when analyzing the movement of crawling cells, for here the con- 
trol of the force-generation machinery is an essential component of the response. 
While amoeboid cells such as Dd extend pseudopods "randomly" in the absence of 
signals, the direction of extension is tightly controlled in the presence of a directed 
external signal, and the direction in which forces are exerted on the substrate is 
controlled via the location of contacts with the substrate. Therefore it is appropri- 
ate to incorporate the force-generation machinery as part of the internal state, and 
as a first step we condense this all into a description of how the force exerted by a 
cell on its surroundings (and vice-versa) depends on the external signal. In reality 
amoeboid cells are also highly deformable, and a complete theoretical treatment 
of taxis at the single cell level has to take this into account. This is currently under 
investigation but will not be pursued here; instead we only describe the motion of 
the centroid of the cell. However, the following framework is sufficiently general 
to allow distributed internal variables within a cell. 

Hereafter we use y as it appears in d!9l > to denote the chemical variables in- 
volved in signal transduction, control of actin polymerization, etc, and we denote 
the force per unit mass on the centroid of a cell by ^(x, v, y). Therefore the inter- 
nal state equations are given by 



dy 

dt 

and the velocity evolves according to 



= G(y,S) (19) 



— =^(x,v,y). (20) 

Here Q : Y x § — > Y is in general a mapping between suitable Banach spaces and 
T : R N x R N x Y -> R N where N = 1, 2, or 3 is the dimension of the physi- 
cal space. This generality is needed because the variable y can include quantities 
that depend on the location in the cell or on the membrane, and which may, for 
example, satisfy a reaction-diffusion equation or another evolution equation. 

The cell is therefore described by the position and velocity of its centroid, 
and the internal state y G Y. In some important cases described later there is a 
projection V : Y — » Z C Y from Y onto a suitable finite-dimensional subspace Z, 
obtained for example by considering the first few modes in a suitable basis for Y, 
such that 

V{Q{y,S)) = G(z,5) and jF(x,v,y) = F(x,v,z), where z = Vy. 

(21) 

Here G(-, S) : Z -> Z and F(-, ■, •) : R N x R N x Z -> R N are mappings between 
finite-dimensional spaces. The first equality defines the function G, though it may 
of course be difficult to find when Q is nonlinear. The function F is explicitly given 
by the second equality when the reduction is possible. 

Given a suitable choice of the projection V, we can reduce the infinite-dimen- 
sional system fl!9l > - ( I20l > to the following set of ordinary differential equations in 
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finite dimensions for the description of individual cells. 

G(z,S) (22) 

F(x,v,z) (23) 

Next, let p(x, v, z, t) be the density of individuals which are at point x, with ve- 
locity v and with the vector of reduced internal variables z; then the transport 
equation dl 31 can be written in the form 

|| + V x -vp + V„-Fp + V 2 -Gp= -A(z)p + J A(z)T(v,v',z)p(x,v',z,i)dv'. 

(24) 

A crucial assumption for using the transport equation formalism is that the pro- 
jection V exists; at present we do not know how to extend this framework to an 
infinite-dimensional manifold. Examples of models for which the projection V 
can be found will be given in the following sections, and in these cases we can use 
(I24t as the starting point for obtaining macroscopic equations. As described ear- 
lier, the right-hand side models the instantaneous changes of direction of motion, 
and in the present context we use this to describe the small fluctuations due to ran- 
dom "errors" in the sensing of the signal and possibly to an intrinsic mechanism 
for random exploration of the local environment. Tranquillo and Lauffenburger 
[50 1 developed a model of amoeboid movement that focuses specifically on the 
stochastic component. 

A natural question is what can be done if a suitable projection V is not eas- 
ily computed, or if the explicit form of G is impossible to obtain because of the 
complexity of the mapping Q. In some cases it may still be possible to describe 
the macroscopic-level dynamics by the evolution of a few slow variables, and by 
using computational equation-free methods which are currently being developed, 
to obtain populational level quantities without explicitly deriving the macroscopic 
equations (see I29I16I15 I and references there), using either the full model of the 
amoeboid cell or the best available reduction of it. 

In the remainder of the paper we give examples of the reduction of d!9i - 
(I20> to the form J22I - \22>\ and the derivation of macroscopic equations via the 
transport equation (I24> . in order to understand how the population-level dynamics 
depends on the characteristics of the individual behavior. We start with a motivat- 
ing example in which we further reduce the system (1221 - d23i by assuming that 
dimZ = 1 and that the function G transduces the signal directly, i.e., z = S. 

2.1. A motivating example 

To illustrate how the effect of acceleration of the cell can enter into the macro- 
scopic equations, we consider the example of motion of a cell in the plane in 
response to a wave of a chemotactic substance. The typical response of Dd or 
leukocytes to a pulse-like wave of the chemoattractant can be divided into several 
phases depending on the position of the cell relative to the wave |21 1. In Figure^ 
we distinguish five different phases - denoted (A) - (E). Before the wave arrives 



dz 
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di 
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direction 
of wave 



Fig. 1. The notation for the different phases of the wave of chemoattractant seen by a cell at 
a fixed spatial position, as a function of time. The horizontal axis is the time and the vertical 
axis is the amplitude of the signal. 

at the cell, there are no directional cues in the environment and the cell extends 
pseudopods in all directions - Phase (E). When the wave arrives the cell experi- 
ences an increasing temporal gradient at all points of its surface and can detect a 
front-to-back spatial gradient over its length (where front denotes the direction of 
the oncoming wave), which causes it to polarize in the direction of the oncoming 
wave. This is Phase (A) in Figure ^ In Phase (B) lateral pseudopod formation is 
suppressed and the cell moves more-or-less directly towards the aggregation cen- 
ter at a speed of 10-20 ^m/min. In natural cAMP waves the cAMP concentration 
at the peak of the wave is high enough that in Phase (C) the cells stop translocat- 
ing and depolarize. In Phase (D) the temporal gradient is negative, although the 
spatial gradient is positive in the outgoing direction, and the cell begins to form 
pseudopods in all directions. This is presumably due to slow adaptation to the de- 
creasing cAMP signal, and as we shall see, if it is too fast the cells may reverse 
direction and follow the outgoing wave. In Phase (E), there is no extracellular sig- 
nal present and there is not net movement of cells. This last phase is not described 
in 1211 but it is of interest to include this to describe the motion in the absence 
of a stimulus. Formal rules used in the context of an individual-based model of 
Dd aggregation show that population-level aspects of chemotaxis such as stream 
formation can be reproduced if the foregoing phases are properly incorporated 1 9 1 . 
How to incorporate these characteristics into a continuum description is the ques- 
tion addressed here. 

The following example is not meant to provide a realistic description of taxis, 
but rather to motivate the analysis done later. In a coarse or high-level description 
of movement in response to signals, information carried by the external signal 
detected by a cell is transduced through the intracellular signaling network, and 
during deterministic turns the velocity of the cell follows the external gradient 
with some delay. We write Newton's law for the motion of the centroid as 



where we assume that the relaxation (or adaptation) time t v (S) is a functional 
of S. Term g(V5) can be interpreted in terms of the force generated from the 



dv _ g(Vff) - v 

dt ~ Tv (s) 



(25) 
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extracellular signal. Typically g(VS) vanishes at zero, is monotone increasing, 
and saturates for large VS. The dependence of t v (S) on S could arise, for instance, 
from different responses of the intracellular dynamics to increasing and decreasing 
signals; or from alterations in the adhesion sites between cell and substrate. In 
earlier work the turning behavior was incorporated via rules 1 9 1, rather than via an 
equation of motion such as d25i . 

To demonstrate that this model can capture some of the salient features of Dd 
aggregation in response to cAMP waves from a pacemaker center, we present the 
results of cell-based numerical simulations that use J25b for the velocity, given a 
suitable choice of t v (S). We consider a two-dimensional disk (corresponding to a 
Petri dish) of radius 5 mm, and we specify a periodic source of cAMP waves at 
the center of the domain. The period of the waves is seven minutes, their speed 
is 400/im/min, and the maximal speed of a cell is about 20 /im per minute, all of 
which are chosen to approximate natural waves in a Dd aggregation field. More 
precisely, we choose g(V5) = sqVS/(c s + \\VS\\) where sq = [20/im/min], and 
c s measures the sensitivity of the signal transduction mechanism. In the numerical 
examples we choose a wave with maximum \\VS\\ equal to 1 mm -1 and c s 
1 'inn] Initially the cells are distributed uniformly and we investigate under 
what conditions the cells aggregate at the source of the waves of chemoattractant 
S. We consider the following two choices for the dependence of t v (S) on the 
external field. 

(1) t v (S) is a constant independent of the signal (cf. Figure[2j. In this case there 
is no aggregation, and in fact, cells move to the boundary of the Petri dish. This is 
not surprising, because cells move in the direction of the increasing gradient of the 
attractant, and cells first move toward the source and then turn around. Although 
the wave is symmetric, the cell movement creates a cellular "Doppler effect" in 
that there is an asymmetry in the time the cell detects the inward-directed gradient 
at the front of the wave versus the time it sees the receding wave. Thus in every 
cycle it moves away from the source longer than it moves toward it, and cells 
eventually accumulate at the boundary. 




Fig. 2. Simulation of 5000 cells that move according to J25II when the relaxation time 
t v (S) is constant. We plot the positions of cells at t — (left) and at t = 4000 min ( right). 
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(2) In this case the relaxation time is specified as a function of the time derivative 
of S at the position of the cell, i. <?., t v (S) = t v (St). t v is chosen so that cells turn 
rapidly when the temporal derivative is positive and slowly when it is negative. In 
our numerical example, we simply put t v = 0.5 min for St > 0, and t„ = 10 
min for St < 0. The results are shown in Figure |5J here one sees that the cells 
aggregate at the source of the waves. 

These cases show that reorientation that is adaptive with respect to the tempo- 
ral gradient of the signal suffices to produce aggregation, as was found earlier in 
formal cell-based rules 1 9 1 and used previously in macroscopic descriptions based 
on the classical chemotaxis equation H45I25I . 




Fig. 3. Simulation of cells which move according to i25ti when the relaxation time T V (S) = 
Tv(St) is chosen so that t v is small (0.5 min) when St is positive and large (5 min) when 
St is not positive. The positions of 5000 cells at different times are plotted. 
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Next we address the derivation of a macroscopic description from the transport 
equation, using the direct effect of the signal on the turning given by l|25}. We 
denote by p(x, v, t) the density of individuals which are at point x G R 2 and have 
velocity v 6 V C IR 2 at time t. Here, V is a bounded, symmetric set which is 
determined by the external signal and by system ( 125 \ . We also assume that there 
is a signal-independent component to the turning for which the kernel T is given 
by T(v, v') = (27Ttio) _1 <5(|v — v'| — Vq), where vq > represents the magnitude 
of the random component of v. The cells add a small random component to their 
velocity at a rate A. Now p(x, v, i) satisfies the transport equation 



dp „ 



p 



(26) 



g(VS) - v 
t v (S) 

Ap(x,v,t) + ^- / <y(|v-v / |-« )p(x,v / ,*)dv'. 



(27) 



2ttw Jv 

We define the macroscopic density n and macroscopic flux j via 

n= p(x,v,t)dv, j= / vp(x,v,i)dv, 

and by integrating J26t over v, and multiplying d26l by v and integrating over v, 
we obtain the following evolution equations for n and j: 

^ + V x -j = (28) 
^ + V x -j-^-ng(VS) = 3 7 -. (29) 

T V (S) T V {S) 

The convective flux — J v ViV^pix, v, t)dv that appears in ( I29l l introduces a 
higher-order moment, and in earlier work on bacterial chemotaxis we could jus- 
tify the closure hypothesis jo- — s 2 nSik/2, where s is the speed of a bacterium 
(cf. 1 18 1). Since the speed is not constant during "runs" in the amoeboid case, we 
must use a different approach here. By constructing the evolution equation for jik 
and assuming that it relaxes rapidly in time, i. e., neglecting its time derivatives, 
and neglecting the third order velocity moments, we find that the 2x2 tensor j 
has components 

r (S)Xv 2 1 

ji*(x,t) = — -nS ik + -(gdk+gkji), fori, fc = 1,2. (30) 

This leads to two possible closures, (i) by keeping only the zero-order moment 
(the term involving n), or (ii) by keeping both the zero-order and first-order contri- 
butions. We use the first of these here and find that the system J28b - ( I29l > becomes 

dn 

7 j- + V x -j = (31) 
at 

<9j (t v {S)\vI \ 1 j 
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In this form we identify the chemotactic velocity and the chemotactic sensitivity 



as 



1 ; 8 (VS) x - 1 8(VS » 



T v (S) ay ' " T V (S) VS 

where the latter only makes sense if we assume that g(VS) = g(VS)WS. If in 
addition g saturates for large arguments, the velocity saturates and the sensitivity 
goes to zero in the presence of large gradients, as one should expect. 

One sees that at this level of closure, the relaxation rate of the flux on the right 
hand side of d32t is signal dependent, but if we were to suppose that r v (S) = tq 
is independent of S then the system (13 II - d32t can be written as the second order 
equation 

d 2 n 1 dn t \v 2 / 1 



8? ' r dt- 4 An-V x .^-ng(V5)j, (33) 

which is the hyperbolic form of the classic chemotaxis equation. However, if t(S) 
is signal dependent the system J3 It — J32l > does not reduce to J33I . and this suggests 
that one cannot expect to obtain the classical form of the chemotaxis equation 
when internal states are taken into account explicitly. On the other hand, as we saw 
in the simulations above, one cannot avoid the "back-of-the wave paradox" without 
a signal-dependent t v 191 121 . Let us note that we can treat similarly a modification 
of i25i where we allow the force to depend directly on the time derivative of the 
signal. This can be done by replacing g(V5) with g(V5, St). 

To illustrate the validity of the macroscopic equation i33\ . let us consider the 
cell-based numerical simulations of Nq individuals whose velocity is governed by 
d25i . We denote the positions of individuals as Xi(t), i = 1, . . . , A*o. Then the 
quantities of interest are the mean position of individuals and the mean square 
deviation, and for the discrete-cell analysis these are defined as 

, N N 

X(^-^x,(t) and a 2 (i) = — £|Mi)-X(t)|| 2 . (34) 
i— l %—i 

By multiplying d33l > by x and integrating over x, then computing the variance, 
much as in Section 1.1, we find that for long times the macroscopic description 
predicts that 

X(i) k, g(VS') t and a 2 (t) rj r 2 Aw 2 t. (35) 

To compare the theoretically-derived results d35i with the cell-based computations, 
we choose tq = 1 min, A = 1 min -1 , vq = 1 ^m/min, iVo = 10 4 cells, g = uld, 
where u> = 20 /im 2 /min, and V5 is given by Q with ||V5||= 1 /itn -1 . We place 
all cells at [0, 0] and set their velocities to initially, compute their subsequent 
motion, and plot the first component of X(t) and a 2 (t) (as given by ( 1341 ) in Figure 
We see that after an initial transient period both quantities grow linearly with 
time, and the slopes are asymptotically equal to the slopes predicted theoretically 
using i BSY 
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(a) (b) 




time [min] time [min] 



Fig. 4. Time evolution of statistics i34l obtained from the cell-based simulation, (a) First 
component of X(t) as a function of time, (b) Mean square displacement a 2 (t) as a function 
of time. 



2.2. The infinite- dimensional model and its finite- dimensional reduction 

As discussed earlier, analysis of E.coli chemotaxis shows that the microscopic 
behavior can be translated into the macroscopic parameters, and it is desirable 
to do the same for amoeboid chemotaxis. However, as noted earlier the internal 
state may now live in a Banach space, and a reduction to finite dimensions is 
necessary. We start with the description of excitation-adaptation dynamics on the 
cellular membrane to model directional sensing and reduce the resulting system. 
For simplicity we suppose that a cell is a disk of radius d. The state of a cell will 
be described by the position x and velocity v of its centroid, and several internal 
variables on the membrane. The membrane of the cell can be described as the set 

M ={d [cos(0), sin(0)] | 9 E [0, 2tt)}. (36) 

The local state at each point of the membrane will be specified by the (infinite- 
dimensional) internal state variable y(6,t) = [yi(9,t),y2(9,t)] T , 6 e [0, 2ir), 
whose evolution is governed by the "excitation-adaptation" cartoon model illl . 
In the formalism of equations d!9l > - i20\ this means that the internal state y(t) : 
[0, 2tt) — ► M 2 can be viewed as an element of the Banach space of 27r-periodic 
vector functions Y and evolves according to 

^(6,t)=S(9,t)T-Ty(9,t) (37) 
for 9 £ [0, 2tt) and t > 0. Here 

S(9,t) = S(x + de(9),t), e(9) = [cos6»,sin6»] T , r = [t" 1 .^ 1 ] 2, , 
and 



The y-variables correspond to those in dl9> . and we project these to finite dimen- 
sions by considering the first Fourier mode of yi and the first two Fourier modes 
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of y 2 . Thus we define the average internal variables as 



1 f 27T 

z(t) = ( Zl (t),z 2 (t)f = — / y(9,t)d0, (38) 



2tt 



o 



<l(t) = ( qi (t),q 2 (t)) T = [ e(9)y 2 (9,t)d6. (39) 
dir Jq 

To derive equations for the reduced finite-dimensional set of internal variables z(t) 
and q(t), we use the approximation 

S(9,t) ~ S(x, t) + d e(0) • VS, (40) 

and consequently, we can write 

dV2 m ,, S(x,t)+de(9)-VS-y 2 (9,t) 

-dt m = 7 a • (41) 

for 9 e [0, 2vr) , t > 0. Multiplying (EJ by 1, cos(0) or sin(#) and integrating the 
resulting equations with respect to 9, we obtain 

dZ2 = 5(x, t) -z 2 
dt r a 
dq = V5(x,t)-q 
dt r a 

Thus Z2 relaxes to the signal 5(x, t) and q relaxes to the directional information 
of S, both with the decay rate r a . To interpret z%, we assume fast excitation (i. e., 
T e = 0). Then using the fact that yi = S — y 2 and integrating J4T1 with respect to 
0, one finds that 

^ = 5 X)i )|v.V%f)^, (44) 
dt at r a 

By integrating this one sees that z\ tracks the Lagrangian derivative of S taken 
along the cell's trajectory, with a memory determined by r a : the smaller r a the 
faster the cell forgets the history of this derivative. Taken together, the four vari- 
ables (zi, z 2 ,qi,q 2 ) contain information about the rate of change of the signal 
along the trajectory (zi), the local value of the signal (z 2 ), and the gradient of the 
signal (q). To this set we will add the polarization axis u = (u\, u 2 ) in the next 
section, and the result will be the smallest set of variables that is able to capture 
the phenomena described earlier. Consequently, the simplest hypothesis is that cel- 
lular motility depends only on these six variables, i. e., the z used in (122b has six 
components. 
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2.3. The motility model 

The next step is to build this model for the internal dynamics into a description 
of cellular movement in order to reproduce some of the experimental behaviors 
observed for eukaryotic chemotaxis. As we saw in Section I2TTI Dd or leukocytes 
respond to the waves of chemoattractant by moving toward the source of waves, 
and the five different stages of the wave with different behavioral responses of 
the cell are schematically shown in Figure ^ These and various other cell types 
also polarize after sufficient exposure to a directional signal 1 26 1 . In order to build 
directional sensing, polarization and response to waves into the model, we dis- 
tinguish three distinct states of cells: (1) polarized cells which are motile (MPC), 
(2) polarized cells which are resting (i. <?., non- motile, denoted RPC), and (3) non- 
polarized cells which are resting (RUC). The signal transduction machinery for 
all types is described by the membrane-based model (I36> - d37l >: the difference 
between the types is in their motility behavior. 

We describe a motile polarized cell (an MPC) by its position x, its velocity v 
and its internal state y, as before. However, instead of defining the force directly 
in terms of the signal, as was done in d25i . we assume that the force is proportional 
to the projected internal variable q (defined by (I39». which tracks the gradient of 
the signal (cf. 1431 ). Thus we write the equations of motion for a cell as 



dx dv 7q — v 

dt dt T4 



(45) 



In a steady gradient of the signal q relaxes to VS on the time-scale r a , and v 
relaxes to 7q on the time-scale 7^; thus the models predicts steady motion in a 
constant gradient. One expects that in general r a < t&. However, as we saw earlier, 
to explain the back-of-the-wave behavior 12 IS the response to the wave must be 
biased toward moving when the signal is increasing in time, as in the front of the 
wave. It is known that Dd cells and leukocytes stop translocating and lose their 
polarity in Phase (C) of a wave (cf. Figure^, which introduces an asymmetry into 
the response, and to capture this we introduce a resting state. A resting cell (either 
an RPC or an RUC) is described by its position x and its internal state y G Y, 
and these cells may also have a polarization axis u = (ui, 1*2). We assume that 
the position of a resting cell is fixed and that the internal state evolves according 

to 07). 

Finally, we must postulate how transitions between the three states depend on 
the signal (cf. Figure[5}. We assume that the motile cells retain their polarity upon 
stopping, and that the transition rate from the motile to the resting polarized state 
depends on Z\ as shown in Figure [5] A moving cell "computes" the directional 
vector q and the average around the perimeter of the internal variable yx, which is 
Z\, according to d43t and (144b . The interpretation of Figure [5] and the justification 
for the postulated dependence of the transition rates between states on Z\ are as 
follows. 



(i) If the Lagrangian derivative of the signal along a cell's trajectory is negative, 
then z\ decreases and the transition rate from the motile state to the resting 
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RUC 




MPC RPC 

Fig. 5. 77ie allowed transitions between cell states. The transition rates depend on the in- 
ternal state z\ as follows: k'21 = Ai — bizi, ki2 = A2 + 62-21, fci3 = A3 + b^zi, /C32 = Ao. 

polarized state, increases. The resting polarized cell adopts a polarization equal to 
the velocity vector before it stops, i. e., u = v after the transition. 

(ii) If a cell is resting and there is an increase in the signal, then z\ increases 
and it is more likely to move. It the cell is unpolarized, than its initial velocity 
(polarization) is zero and the time reflects the time delay needed for polarization 
of the cell. If the cell was already polarized than its initial velocity is equal to the 
polarization vector, i. <?.,. we set v = u, and the time delay reflects the relaxation 
time for turning, if it is necessary. 

(iii) A resting polarized cell looses its polarity at a rate Ao, and thus polarized cells 
which do not receive a stimulus for a long time lose their polarity. 

Next we demonstrate that the model can successfully solve the back-of-the-wave 
problem 1 121211 . i. e., cells will aggregate at the source of the attractant waves. 

2.4. Aggregation when resting states are incorporated 

The internal dynamics model y written in terms of dl 91 is given by d36l - ( I37> . and 
every cell is described by its position x£l 2 , its velocity v G R 2 , its polarization 
axis u £ R 2 and its internal state function y 6 Y. To compute z\ and q, the radius 
of a cell is set to d — 7.5 /im, and we discretize the cell boundary (13 6> using m 
meshpoints, 

9j = — -. for j = 1, 2, . . . , m — 1, to. 

to 

Then the state of each cell is described by an (to + 4) -dimensional vector 

(^,y(0 1 ),y(e 2 ),...,y(0 m )). (46) 

Here, v denotes the velocity for an MPC and the polarization axis for an RPC. 
We can simply set this equal to for RUCs, which are in fact described by to + 
2 variables. The internal state variables y(6j) evolve according to m equations 
of the form (1371 . which are uncoupled because there is no transport along the 
membrane. The evolution of x and v is described in Section 1231 At each time 
step we use the y(9j ) to numerically approximate integrals J38I - (1391 and thereby 
compute z\, which is needed for the computation of the transition rates between 
different states as shown in Figure|5] and q, which is necessary for the integration 
of ( 145 \ . Throughout we use to = 50, and therefore every cell is described by a 
54— dimensional state vector ( 1461 . 
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As was done for Figures |2] and [3] we consider a two-dimensional disk of ra- 
dius 5 mm, and we specify a periodic source of cAMP waves at the center of the 
domain. The period of the waves is seven minutes, their speed is 400 /im/min, and 
the waves are scaled so that the maximum ||V5|| is 1 mm -1 . Initially the cells are 
distributed uniformly. We use the following base transition rates and sensitivities 
in the transition rates fey given in Figure|5J 

Ai = A2 = A3 = 1 min -1 , Ao = 0.2 min -1 and b± = 62 = &3 = b (47) 

Later the parameter b will be varied (cf. Figures[6]and[7ll. The time constants are 
chosen as 

r e = (fast excitation), r Q = 0.5 min, = 2 min. (48) 

The two parameters which have yet to be specified are 7 in equation J45> and b. 
The parameter 7 simply rescales the speed of cells. We know from experiments 
that the maximal speed of a cell is about 20 /im per minute, which can be used to 
fit the value of the parameter 7. We found that for 7 = 0.08 mm 2 /min, the average 
speed of cells on the steepest part of the wave front is between 10 /im per minute 
and 20 /im per minute in all simulations. Hence, we used 7 = 0.08 mm 2 /min to 
compute the plots shown in Figures[6]and0 

The parameter b specifies how strongly the turning rates depend on z\ and we 
tested three possibilities b = 0, b = 1 min -1 and 6 = 2 min -1 . If b = the transi- 
tion rates kij are independent of z\, and the time evolution of the cell positions is 
shown in Figure[6] We see that in this case there is no aggregation, which is similar 
to what was shown earlier in Figure[2]where we considered the model without the 
internal dynamics and with a constant relaxation time. The computational results 




Fig. 6. The cell distribution as a function of time for b = 0. Positions of 5000 cells at times 
min and 1000 min for periodic waves of chemoattractant. 

for 6 = 2 min -1 are shown in Figure |7] where the aggregation time is compara- 
ble to the eight hours observed experimentally. The results for 6=1 min -1 are 
similar - the only difference is that the aggregation is slower (results not shown). 
Using (I47> . we see that the transition rates (which depend on z\) can be expressed 
in units of min -1 as fc 12 = fc 13 = 1 + bz\ and fc 21 = 1 — bz\. Since bz\ is 
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approximately in range [—0.35, 0.35] min 1 for b = 1 min 1 and in the interval 
[—0.7, 0.7] min -1 for b = 2 min -1 , it implies that the turning rates are in the inter- 
val [1 - 0.35, 1 + 0.35] min -1 for b = 1 min -1 and in the interval [1 - 0.7, 1 + 0.7] 
min -1 for 6 = 2 min -1 . As will be seen in Section|5J the moment approach used 
there is justified when bz\ is small, i. e., for small bias of the turning rates. The 
error may increase significantly for large b because the higher order moments may 
not be negligible. 




Fig. 7. The cell distribution as a function of time for b = 2 min~ . Positions o/5000 cells 
at times min, 100 min, 200 min, 300 min, 500 min and 1000 min for periodic waves of 
chemoattractan t. 

In these figures, as in the preceding ones related to aggregation, there is no 
stream formation such as is observed during aggregation of Dd. Here cells always 
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move radially inward toward the source because the waves are imposed and are ax- 
isymmetric. In the presence of signal relay, as in Dd, it was shown earlier 1 32 1 that 
signal relay combined with a random initial distribution of cells plays an essential 
part in stream formation. 



3. Transport equations 



Next we show that the microscopic model for signal detection, transduction and 
movement can be embedded in a system of transport equations and thence into 
a system of moment equations for macroscopic quantities. To that end, note that 
every cell with the same x, v, Zi,q and same polarization state will follow the 
same rules for movement, so it is natural to introduce density functions p 1 , p 2 , p 3 
as follows: 



• p 1 (x, v, z\ , q) is the density of moving cells at position x with velocity v and 
internal moments z\, q; 

• p 2 (x, u, z\ , q) is the density of resting polarized cells at position x with po- 
larization axis u and internal moments Z\, q. To simplify the form of resulting 
transport equations, we denote the polarization axis as u = v in what follows; 

• p 3 (x, z\ , q) is the density of resting unpolarized cells at position x and with 
internal moments z\, q. 

Hereafter we assume excitation is fast, i. e.,T e = 0, and we use the approximation 
given at (1401 for the signal. The evolution of internal variables q, z\ is therefore 
given by (I43> - J44i . In order to simplify the following equations for p 1 , p 2 and 
p 3 , we define an operator L by 



Cr 



Or 

at 



(VS - q) r 



d 


\(--~)r] 


dzi 


\0t Ta) 



(49) 



then the transport equations for p 1 , p 2 and p 3 are 



Cp 1 = -V x • vp 1 - v v • 



Td 



( 7 q- v)^ 1 



dzi 



[(WS)p 1 } 



-(Ai - b lZl )p l + (A 2 + b 2 z 1 )p 2 + (5 V (A 3 + b 3Zl )p 3 , 
Cp 2 = (Ai - bizijp 1 - (A 2 + b 2 z 1 )p 2 - A p 2 , 

£p s = \ [ p 2 dv-(X 3 + b 3Zl )p 3 . 
Jv 



(50) 
(51) 

(52) 
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where S v is the Dirac function. Next we define the macroscopic densities of parti- 
cles in different states as follows: 

n 1 (x,t)= / / p 1 (x,v,zi,q)dvd2;idq, (53) 
Jv Jr 3 

/ / p 2 (x,v,zi,q)dvdzidq, (54) 
Jv Jr 3 

/ p 3 (x, zi,q)dzidq, (55) 

Jr 3 



n 6 (x,t) 

n(x, t) = n 1 (x, t) + n 2 (x, t) + n 3 (x, t) , (56) 

where n(x, t) is the total density of cells. Here and hereafter the superscript i 
denotes a quantity associated with the i th species, for i = 1,2, 3. If an evolution 
equation in n(x, t) alone could be found the problem would be reduced to the 
classical case. However we will see that this is not possible in general. First, we 
define some additional moments that arise in the usual manner from J50i - d52i 
during derivation of moment equations. More precisely, we derive the evolution 
equations for J53I - (I56> and for the following moments 

ii(x,t)=/ / w fc p l (x,v,zi,q)dvdzidq, £ = 1, 2; fc = 1, 2; (57) 



/ / v k p l (x,v, zi,q)dvdzidq, 

Jv Jm. 3 

n LCM)=/ / %P*(x,v,^i,q)dvd2;idq 
Jv Jm. 3 



i = 1,2; k = 1,2; (58) 



n*( X) t) 
<(x,t) 
n s z (x,t) 



/ / zip l (x, v,zi,q)dvdzidq £ = 1, 2; (59) 

Jv Jr 3 

/ q fc p 3 (x, zi,q)dzidq, (60) 

/ zip 3 (x, zi,q)dzidq. (61) 

JR 3 



Multiplying (I50i - d52i by 1, v\, V2, qi, q2 or z and integrating with respect to v, q 
and zx, we obtain the evolution equations for moments (1531 - i61\ . This system of 
partial differential equations is not closed - it contains some higher order moments 
of the following form 

mUtotoMte ( x ' f ) =/ / Wi 1 ^2 2 ?i fe3 ?2 4 ^ 5 f i (x,v,2 1 ,q)dvdqdzi, 
JvJm 3 

ml uk2 u 3 (x,t) = / < 7 * 1 ^ 2 z 1 fc V(x,v,z 1 ,q)dqd0i, 
Jr 3 

where i = 1, 2, fc Q , a = 1, 2, 3, 4, 5, are nonnegative integer, and the superscript 
fc a on terms in the integral denotes the fc Q -th power of the corresponding vari- 
able. The simplest way to close the moment equations is by setting to zero all 
higher-order moments which do not appear in J53i - (16 11 . More precisely, we use 
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the following closure assumption: = m\ 1 = m\ = toJ 2000 = 
m i, 0,1, 0,0 = m o, 1.1,0.0 = TO i,o, 0,1.0 = m o. 1,0, 1,0 = TO i,o.o,o.i = m o. 1,0. 0,1 = 

\ 1 x 2 2 '' 2' 

m 0, 0,1:0,1 — m 0, 0,0, 1,1 — m o, 0,0, 0,2 — TO 1,0,0,0,1 — m 0, 1,0, 0,1 — m 0, 0,1, 0,1 — 

TO o,o,o,i,i = TO o,o,o,o,2 — TO i,o,i — TO o,i,i — TO o,o,2 = 0- This closure assump- 
tion can be justified for shallow gradients of the signal 1171 . 

Under this assumption, we multiply J50i - d52l by 1, V\, v 2 , <h, 92 or z, we 
integrate with respect to v, q and Z\, and we discard the higher-order moments; 
the result is the following closed system of 16 macroscopic equations. 

°" 1 + ^ + ^ = -Am 1 + A 2 n 2 + A 3 n 3 + b x n\ + b 2 n\ + b 3 n 3 z , (62) 



dt dx\ dx 2 
dn 2 
~dt 
dn 3 
~dt 



= Xin 1 - (A 2 + \ )n 2 - b x n\ - b 2 n 2 z , (63) 
= X a n 2 - X 3 n 3 - b 3 n 3 z , (64) 



Q-f- Trf \-!'"qk J vk\ " LJ v k ' " Z -Jv k 



inl - jl] = -Aid + X 2 j 2 , k = 1, 2, (65) 



dj 



■2 



dt 



Ai4 -(A 2 + A )^ fc , k = l,2, (66) 



dn 1 

— -—- 1 i-~ 1 --\in 1 qk +\ 2 n 2 qk +\ 3 n 3 qk , k = l, 2,(67) 



dt 
dn 3 



dt 



1 

T a 


dS , 


1 


Ik 


1 


dS , 
^— n ^ 


1 


n lk 


r a 




1 


as 
^— n ^ 

dxk 


1 


< 


T a 





Am 1 -(A2 + A0X, fc=l,2, (68) 



' ' ■ • - A r4-A 3 r4, fc = l,2, (69) 



"af " W n + ^ " E dirA = - Xinl + X2nl + Xsnl (70) 

i—1 

= Xml - (A 2 + A )n 2 , (71) 



9n 2 dS 9 1 9 



— n 3 + — n\ = \ n 2 z -X 3 n 3 z (72) 

dt dt T a 

Note that the sum of equations j62t - j64t is the standard continuity equation for 
n, /. e., 

^^ + ^=0. (73) 



<9i 9#i dx 2 
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The system i62\ - 1721 can be written more compactly by defining the following 
vectors and matrices. 



n = (n 1 ,n 2 ,n 3 ) T 



n z = (nl,n 2 z ,nl) T 



n q = (n 9l , n g2 ) = (n 1 , n 2 n 3 , n 1 , n 2 , n 3 v ' 



qi ' qi qi ' 92 » 92 ' 92 / 



j (jt>l J Jt>2 ) Oil! i Jux ' Jt) 2 ' 3v 2 ) 



V 



_a_ _d_ 



A = 



B 



-Ax 


A 2 




A3 " 




Ai 


-(A 2 + A ) 





A 





A 




-A3 




61 


h h 






"1 


-61 


-62 






J = 








-63 










— Ai A2 
Ai -(A2 + A0) 



Ji 



1 




(74) 



(75) 



We further define the tensor product of an s% x S2 matrix X = {xikYi^kJ^i w ' trl an 
S3 x S4 matrix Y = {yikYi fc=i t0 be the (S1S3) x (S2S4) matrix 



X (g) Y = 



cci ; iY xi j2 Y . . . xi, S2 Y 



Then J62t - J72i can be rewritten in the form 



9n 



[V J) j = An + Bn z , 

dn z as ( 1 



V J S8Jj, 



(76) 
(77) 



9n 9fc _ (VS) fc 



n+ [A 

T„ \ Tn 



1 



I3 I n„ 



3T - T/^ 



A± Ji )j Vk 



for k = 1,2, (78) 
for fc = 1,2. (79) 



wherein Ij. is fc x fc identity matrix. This system can in turn be written more com- 
pactly as the system 

,>V DU = A(x,t)U (80) 



Of 



wherein 



U = 







"0 




















n 2 


D = 














n q 
















V j J 

















n = v T ®j, (8i) 
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and 



A 



8S T 



-!l3 

Ta 





Td 



A -I 3 

h ® J T 



(vs) rj 



Jl 



(82) 

We should note that discarding higher-order moments can be justified for the small 
signal gradient case in which the second moments of internal variables are suffi- 
ciently small compared to lower- order moments 1171 . It is important to note that 
the second order velocity moments m\ 1 000 , m\ an( ^ m o 2000 were a ls° 
set to zero because we do not have an obvious moment closure for them similar to 
what was used in the bacterial case [ 18 1, where the Cattaneo approximations could 
be used. 

To obtain a better approximation of these moments we can follow the reasoning 
that lead to the closure d30i earlier. To illustrate this, let us modify the taxis model 
by adding a random component to the cellular movement, namely we change the 
transport equation J50i to the following equation forp 1 : 



Cp 1 = -V x • vp 1 - V v 



Td 



( 7 q- v)p x 



dzi 



[(WS)p 1 } 



-Xp 1 + X / T(v,vV(x,v / ,t)dv / -(Ai-Mi)p 1 
Jv 



(83) 



+ (A 2 + b 2Zl )p 2 + <5 V (A 3 + b 3Zl )p 3 



where the turning kernel is given by T(v, v') = (2irvo)~ 1 5(\v— v'| — vq) similarly 
as before, and vq and A are positive constants. The rationale behind this kernel is 
to incorporate some noise into the system, and vq specifies the "strength" of this 
noise. If we follow the previous procedure, we would obtain the same system of 
equations d80i . which are independent of Vq- This is not suprising - we already 
saw in Section l2~T1 that equation (1331 contains the diffusion term if the apropriate 
closure assumption is derived from d30i and used for the second-order velocity 
moments. Similarly as in d30i . we can multiply the transport equations (I84> and 
(15 1 1 by VkVi, k, I = 1, 2, and neglect time derivatives of the convective fluxes, 
third-order velocity moments and mixed velocity-internal moments to obtain 



in.-, 



AT d (Ao + X 2 )vl 
4(A + A 2 ) + 2r d A A 



-Tl^X, t), 



(84) 



rn 



1,1,0,0,0 



0. 
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Using the moment closure J84l > for the convective momentum flux of p 1 , we obtain 
the following system of moment equations (compare with J80». 



dU 

— +DU = A(x, t)U 



(85) 



Here A(x, t) is given by d82i and U = (n, n 2 , n q , j) T is as in 18 H . but D is now 
given by 



D 



o o o n 
oooo 



a.n T 



where 



a 



Ar d (A + \ 2 )vq 
4(A + A 2 ) + 2r d A Ai ' 



(86) 



3.1. Analysis of the statistics of motion 

To illustrate the validity of reducing the transport equation to the system of hy- 
perbolic equations, we compute the dependence of the mean speed of cells on the 
strength of the underlying signal. 

Multiplying equation (I73> by x and integrating with respect of x, we obtain 
the equation for the mean speed v a „ (t) of the cellular population in the following 
form 



d_ 

dt 



1 



"0 JE 2 



xn(x, i)dx 



no 



where 



n = 



i(x, i)dx and j 1 = / j 1 (x, t)dx. 



(87) 



(88) 



Here, uq denotes the total number of cells in the system and j 1 is the spatial average 
of the flux j 1 = , j% ]. Consequently, to estimate the average speed of cells at 
a given time we have to estimate j 1 /no. 

To do this we use a one-parameter family of time-independent linear distribu- 
tions of extracellular signal S defined by (0, parametrized by |j VS ||. Then the 
matrix A(x, t) = A(|| V5 ||) is independent of x and t, and we can integrate 
equation d80t with respect to x to obtain 



dU 



A(||V5||)U where 



U = 



U(x,<)dx. 



(89) 



Solving system (1891 for U, we can estimate the value of mean speed of the cells 

as 



v Q „(i) = — 
n 



U 



13 



Ui + U 2 + U 3 

and we see that v av (t) will asymptotically approach the velocity given by 

V>13 



(90) 
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signal size ]|VS|] [mm ] time [min] 



Fig. 8. (a) Comparison o/v^ computed by \90\ (solid curve) with results obtained by 
stochastic simulations (circles) for different values of \\\/ S\\. (b) Average position of indi- 
viduals (given by stochastic simulation) as a function of time for ||V 'S\\ = 0.2 mm -1 . 



where if> is a solution of 

A(||VS||)V = 0. 

Using parameter values (147 > - (I48> with 6=1 min -1 and 7 = 0.08 mm 2 /min, 
we can compute the asymptotic average speed by ( I90i for different values of 
||V/f>||. The solid curve in Figure|FJa) shows vj^ as a function of ||V5||. 

The theoretical result d90t can be verified by stochastic simulations, and to 
that end we consider an ensemble of 500 cells. We discretize the boundary d36l 
of each cell using to = 50 mesh points. Hence, the state of each cell is described 
by 54-dimensional vector J46i similarly as in Section 1241 The internal dynamics 
model y written in terms of dl9l is given by d36i - d37i . The radius of a cell is set 
to d = 7.5 //m and we use parameter values (I47> - d48i with 6=1 min -1 and 
7 = 0.08 mm 2 /min. The initial conditions are the same for all cells: namely all 
cells begin at position x = 0, their initial velocities satisfy v = 0, their internal 
variables are equal to around the entire membrane, and cells are initially unpo- 
larized. The average position of cells as a function of time is given in Figure |8jb) 
for ||V5||= 0.2 mm -1 . Since cells are initially unpolarized and resting, the initial 
cellular flux is zero. If we wait for a sufficiently long time, the average speed of the 
cells relaxes to a constant, and when we estimate this from long-time simulations, 
we obtain the values which are shown as circles in Figure|8la). Comparing data in 
Figure|8la), we see that the theoretical result d90i gives a very good approximation 
of the mean asymptotic speed estimated from simulations. This demonstrates the 
fact that one can extract population level information from the moment equations 
derived earlier. 

Finally, we note that the previous analysis can be repeated for d85l . The differ- 
ence between d80i and d85i is the additional noise in the latter, which leads to the 
system (I85> . However, this noise will only influence the diffusion constant and the 
average speed of the population will be unchanged. 

3.2. Further reduction of moment equations 

One can further reduce the size of the system of moment equations d80i by sup- 
posing that the internal dynamics evolve much faster than changes in cytoskeleton 



30 



R. Erban and H. G. Othmer 



and movement. Considering the excitation time r e = and that the adaptation 
time r a -c Td, we can assume the quasi-equilibrium in the equations for n q and 
n 2 in J80l l. i. e., 



T a (I 3 - TaA)- 1 



||„+{(V5f ®J}j 



(91) 



and 



n q = [I 2 ® (I 3 - r a A)] l [VS®n] (92) 

Substituting formulas ( 19 U - ( I92t into ( I80> . we can formally derive the reduced 
system of 7 moment equations for n and j only. These equations are derived under 
the assumption that r a is small. Passing formally to the limit r a — > in ( I91> - 
(I92> . we obtain n z — * and n q — > VS 1 (g> n. As one would expect in light of the 
discussion following (I45> . the reduced equations predict movement up a steady 
gradient, but not in a periodic wave for r a = 0. 

Another approach to eliminate the internal dynamics is to assume the quasi- 
equilibrium assumption directly in J43i - (I44i . i.e. 

dS 

q(x,t)»VS(x,i) and Zl (x, i) fa T Q — (x, t) + r Q v • VS(x, i). 

or 

Denoting p 1 (x,v) the density of motile cells at position x € K 2 with velocity 
v e C R 2 , p 2 (x, v) the density of resting polarized cells at position x with 
polarization axis v and 7i 3 (x) the density of resting unpolarized cells, we can 
write transport equation for p 1 , p 2 and n 3 . Again 7 moment equations for n and 
j can be derived. Such equations were derived and analysed for one-dimensional 
case in 1131 . It was shown that in some parameter regimes, the reduced system 
approximates the simulation with a reasonable precision. See 1 13 1 for details. The 
precision of the approximation depends on the parameter values chosen. 

Finally one can ask whether the method developed in 1171181 . for reducing 
a hyperbolic system similar to d80i to the classical chemotaxis equations, can be 
applied here as well. In the context of chemotaxis based on a "run-and-tumble" 
strategy we were able to analytically compute the eigenvalues and eigenvectors 
of A(x, t), and it was shown that they are independent of the chemotactic signal. 
By exploiting the facts that the spectral gaps are signal-independent, and that rea- 
sonably simple formulas for eigenvectors are available, we reduced the hyperbolic 
system of moment equations to the classical chemotaxis description in the bacte- 
rial case 1 17 18 1. In the case of system i80i . the resulting slow eigenvalues are, for 
general values of parameters, very complicated functions of the parameters, and in 
particular they depend on the signal. To illustrate this, we consider the linear sig- 
nal distribution (0 that leads to ( I89> . we consider parameters from Section lXTl and 
plot the real parts of the eigenvalues of A^ViS]]) as functions of \\VS\\ in Figure 
[9] One sees there that the eigenvalues vary significantly with the signal strength. 

Further analysis is needed to better define the conditions under which the hy- 
perbolic system can be reduced further. In that vein, we note that ignoring the term 
5Zj = i(VS')ij^. in equation J70i leads to the matrix A(|| V5||) which has signal- 
independent eigenvalues. Thus the possibility of further reduction clearly depends 
on the closure assumptions. 
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(a) (b) 




signal size ]|VS|] [mm 1 ] signal size ]|VS|] [mm 1 ] 

Fig. 9. The real parts of the eigenvalues of A(|( VS\\) for (a) ||VS|[e [0, 1] mm -1 ; (b) 
||VS[[€ [0, 10] mm" 1 . 

4. Discussion and conclusions 

The goal in this paper was to derive macroscopic equations for the collective be- 
havior of amoeboid cells based on models for individual-level behavior. In pre- 
vious work we developed a moment closure approach to the transport equation 
for a velocity jump process that describes cell motion for cells that use a "run- 
and-tumble" strategy 1 17 18 1. Here we demonstrated that this approach can also 
be applied to the more complex processes involved in the movement of crawling 
cells, and showed that one can predict important macroscopic characteristics from 
knowledge of the individual-level properties of these cells. We focused on chemo- 
taxis in Dd as the model system because much is known about this system, but 
the general approach can be applied to any type of extracellular signal, including 
those that arise from all receptor-based interactions of a cell with its environment. 
Here we summarize the approach and discuss its advantages and limitations. 

In Section|2]we introduced the general model d!9l > - J20i for the behavior of an 
individual eukaryotic cell. Since this model is often infinite-dimensional, we have 
to first reduce it to the finite-dimensional form J22I - d23l >. If such a reduction is 
possible, we can apply the transport equation framework J24i for the reduced set 
of variables. In this case we can derive the appropriate moment equations J80l > and 
use them to study the macroscopic collective properties of cells, as we illustrated in 
Section im where we studied the dependence of the average speed of the cellular 
population on the strength of the extracellular signal. 

Therefore the crucial assumption for a model which can be treated in the frame- 
work developed here is that the projection V from J2 It exists and the equations 
J22I - d23i can be easily written. If this is not the case, it may be still possible to 
reduce the individual-level dynamics to a low-dimensional description of an indi- 
vidual cell. The behavior of these coarse (intracellular) observables (on the level of 
a cell) can be studied by computational equation-free methods which are currently 
being developed H29I15I . The similar computational methods can be then used to 
study the population-level properties of the amoeboid cells using either the full 
model of an individual cell or the best available reduction of it 1 291161 . 

The moment closure reduction of the transport equation used here can be jus- 
tified for small signal gradients 1 17 1, but in the case of large signal gradients, the 
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higher order moments may not be negligible and cannot be discarded. Similar mo- 
ment methods can be used for any model assuming that the internal dynamics is 
close to its quasi-equilibrium. If the original internal dynamics model is nonlinear, 
it can be linearized around its equilibrium value for small signal gradients and the 
moment approach can be applied. Hence, the reader should view our linear model 
as an example of the linearization of more complex nonlinear problems. Of course, 
the linear model is clearly not valid for large signal gradients, but this is not the 
parameter regime studied in this paper. However, even some strongly-nonlinear 
models for internal dynamics produce simple input-output behavior that can be 
captured by a linear model with possibly signal-dependent parameters [18|, and 
thus the results presented here may have broader applicability than the deriviation 
would suggest if applied strictly. 

In the case of a constant external signal gradient we were able to derive explicit 
expressions for various statistics of the motion from the hyperbolic system derived 
from the transport equation. We also discussed the predicted behavior of models 
for experimental conditions such as spatio-temporal waves of chemoattractant. It 
is known that eukaryotic cells such as Dd or leukocytes aggregate at the source of 
the waves, and the models studied here include the processes, such as adaptation, 
that are necessary to reproduce this behavior. 

The models described here are all based on deterministic extracellular signals, 
although a small random component was added to the choice of direction. The 
effects of stochastic fluctuations in signal detection and processing on movement 
were introduced in II II . where stochastic differential equations are postulated to 
model cell movement on the time scale of the molecular processes that govern 
locomotion. Some concrete estimates of the probable role of noise in the signal 
seen be a Dd cell were made in 1391 . However a much more detailed analysis of 
stochastic effects in all components of the movement response is needed. 
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